
function [dfdt] = waves(t,y,L,C,E)
global Nx Ny;
dfdt(1:2*Nx*Ny,1) = -L*y(1:2*Nx*Ny,1) + C * y(2*Nx*Ny+1:3*Nx*Ny,1);
 G = full(C);
% 
 F = C(1:Nx*Ny,1:Nx*Ny);
 F = sparse(F);


dfdt(2*Nx*Ny+1:3*Nx*Ny,1) = E*y(2*Nx*Ny+1:3*Nx*Ny,1) + F*y(Nx*Ny+1:2*Nx*Ny,1);